function F_star = DCSX_mktwide(df_sim,coefs_sim,dfull_sim, model, nmarket)

brndfe = df_sim.brndid == unique(df_sim.brndid)';
montfe = df_sim.montid == unique(df_sim.montid)';
prodfe = df_sim.prodid == unique(df_sim.prodid)';

% Random sample X of markets
unique_cdid = unique(df_sim.cdid_sim);
sampled_cdid = (1:nmarket)';
% sampled_cdid = randsample(unique_cdid, nmarket);

sel = ismember(df_sim.cdid_sim, sampled_cdid);
df_sim = df_sim(sel, :);
dfull_sim = dfull_sim(sel, :);

% data parameters
fesd     = [prodfe(sel, :) montfe(sel, :)];
x        = [df_sim.p_jt fesd];
xnonlin  = [df_sim.p_jt ones(size(df_sim.p_jt)) df_sim.calor];
xi       = [df_sim.xi];

% model parameters
if model == "L_brandfe_montfe" || model == "RCL"
    fes_M       = [brndfe(sel, :) montfe(sel,:)];
elseif model == "L_brandfe"
    fes_M       = brndfe(sel, :);
end
x_M         = [df_sim.p_jt fes_M];
if model == "RCL"
    param_M0    = zeros(size(x_M,2)+2,1);
else
    param_M0    = zeros(size(x_M,2),1);
end

%create df_sim_wide
x_mrkwide = zeros(size(prodfe, 2), size(x, 2), nmarket);
xi_mrkwide = zeros(size(prodfe, 2), size(xi, 2), nmarket);
xnonlin_mrkwide = zeros(size(prodfe, 2), size(xnonlin, 2), nmarket);
dfull_sim_mrkwide = zeros(size(prodfe, 2), size(dfull_sim, 2), nmarket);
in_samp_mrkwide = zeros(size(prodfe, 2), size(dfull_sim, 2), nmarket);

for j = 1:size(prodfe, 2)
for m = 1:nmarket
    sel = (df_sim.cdid_sim == m) & (df_sim.prodid == j);
        if any(sel)  % Check if there are any matching rows
            x_mrkwide(j, 1:size(x(sel,:), 2), m) = x(sel, :);
            xi_mrkwide(j,1:size(xi(sel,:),2),m) = xi(sel,:);
            xnonlin_mrkwide(j,1:size(xnonlin(sel,:),2),m) = xnonlin(sel,:);
            dfull_sim_mrkwide(j,1:size(dfull_sim(sel,:),2),m) = dfull_sim(sel,:);
            in_samp_mrkwide(j,1:size(dfull_sim(sel,:),2),m) = 1;            
        end
end
end

f_obj = @(param_val) ObjFun(param_val,df_sim,coefs_sim,x_mrkwide,xnonlin_mrkwide,xi_mrkwide,dfull_sim_mrkwide, in_samp_mrkwide, nmarket, model, brndfe, montfe);

   options = optimoptions('fminunc',     ...
        'Display', 'iter-detailed',       ...
        'MaxIterations',300,              ...
        'Algorithm', 'quasi-newton',    ...
        'StepTolerance', 1e-10,           ...
        'SpecifyObjectiveGradient', true, ...
        'OptimalityTolerance',1e-6,       ...
        'checkGradients',false            ...
    );

% Perform optimization
[param_star] = fminunc(f_obj, param_M0, options);

% Return optimized value
[F_star, dF_star] = ObjFun(param_star,df_sim,coefs_sim,x_mrkwide,xnonlin_mrkwide,xi_mrkwide,dfull_sim_mrkwide, in_samp_mrkwide, nmarket, model, brndfe, montfe);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [F,dF] = ObjFun(param_M,df_sim,coefs_sim,x_mrkwide,xnonlin_mrkwide,xi_mrkwide,dfull_sim_mrkwide, in_samp_mrkwide, nmarket, model, brndfe, montfe)

theta    = [coefs_sim.alpha; coefs_sim.sigma_D];
Pi       = coefs_sim.Pi;
rho      = coefs_sim.rho;

J = max(df_sim.prodid);
K = length(x_mrkwide(1,2:end,1));

% Create theta_M, which has dimensions theta
brand_FE = repelem(param_M(2:size(brndfe,2)+1), 3);
theta_M = [param_M(1); brand_FE];
if model == "L_brandfe_montfe" || model == "RCL"
    theta_M = [theta_M; param_M(2+size(brndfe,2):1+size(brndfe,2)+size(montfe,2))];
else 
    theta_M = [theta_M; zeros(size(montfe,2),1)];
end
% random coefficients for model
if model == "RCL"
    Pi_M = [0; param_M(end-1:end)];
    nparam = size(theta_M,1) + 2;
else
    Pi_M = [0;0;0];
    nparam = size(theta_M,1);
end
rho_M = 0;

% Preallocate memory for temporary variables
Delta_Y         = zeros(J,J,nmarket,K);
Delta_Y_M       = zeros(J,J,nmarket,K);
GradDelta_Y_M   = zeros(J,J,nmarket,K,nparam);

% Calculate data shares for data
[share_main, ~]      = f_share_mktwide( ...
    theta,              ...
    Pi,                 ...
    rho,                ...
    x_mrkwide,          ...
    xnonlin_mrkwide,     ...
    xi_mrkwide,  ...
    dfull_sim_mrkwide,   ...
    in_samp_mrkwide      ...
);

% Calculate xi_M for logit
if model == "L_brandfe_montfe" || model == "L_brandfe"
    xi_M = log(share_main) - log((1-sum(share_main,1))) - pagemtimes(x_mrkwide, reshape(theta_M, [], 1, 1));
    % make sure not inf
    inf_indices = isinf(xi_M);
    xi_M(inf_indices) = 0;
% Calculate xi_M for RCL
else 
    [Delta_RCL, dDelta_RCL] = rcnl_delta( ...
        share_main,         ...
        Pi_M,            ...
        rho_M,           ...
        xnonlin_mrkwide,          ...
        dfull_sim_mrkwide          ...
    );

    xi_M = Delta_RCL - pagemtimes(x_mrkwide, reshape(theta_M, [], 1, 1));
end

% Calculate data shares for data
[share_main_M, share_main_ijt_M]      = f_share_mktwide( ...
    theta_M,              ...
    Pi_M,                 ...
    rho_M,                ...
    x_mrkwide,          ...
    xnonlin_mrkwide,     ...
    xi_M,  ...
    dfull_sim_mrkwide,   ...
    in_samp_mrkwide      ...
);

for k = 1:K
for j = 1:J
x1 = x_mrkwide;
x0 = x_mrkwide;

x1(j,k+1, :) = 1*in_samp_mrkwide(j,1,:);
x0(j,k+1, :) = 0*in_samp_mrkwide(j,1,:);

% Calculate shares for zero
if isequal(x0, x_mrkwide)
    share_0 = share_main;
    share_0_M = share_main_M;
    share_main_0_ijt_M = share_main_ijt_M;
else
    [share_0, ~]      = f_share_mktwide( ...
        theta,              ...
        Pi,                 ...
        rho,                ...
        x0,          ...
        xnonlin_mrkwide,     ...
        xi_mrkwide,  ...
        dfull_sim_mrkwide,   ...
        in_samp_mrkwide      ...
        );

    [share_0_M, share_main_0_ijt_M]    = f_share_mktwide( ...
        theta_M,            ...
        Pi_M,                 ...
        rho_M,                ...
        x0,         ...
        xnonlin_mrkwide,     ...
        xi_M,  ...
        dfull_sim_mrkwide,   ...
        in_samp_mrkwide      ...
    );    
end

if isequal(x1, x_mrkwide)
    share_1 = share_main;
    share_1_M = share_main_M;
    share_main_1_ijt_M = share_main_ijt_M;    
else
    [share_1, ~]      = f_share_mktwide( ...
        theta,              ...
        Pi,                 ...
        rho,                ...
        x1,          ...
        xnonlin_mrkwide,     ...
        xi_mrkwide,  ...
        dfull_sim_mrkwide,   ...
        in_samp_mrkwide      ...
    );

    [share_1_M, share_main_1_ijt_M]    = f_share_mktwide( ...
        theta_M,            ...
        Pi_M,                 ...
        rho_M,                ...
        x1,         ...
        xnonlin_mrkwide,     ...
        xi_M,  ...
        dfull_sim_mrkwide,   ...
        in_samp_mrkwide      ...
    );
end

Delta_Y(:,j,:,k) = share_1-share_0;
Delta_Y_M(:,j,:,k) = share_1_M-share_0_M;

% Calculate gradient
for p = 1:length(theta_M)

    % calculate sum for M1
    sum_product_1_M = sum(share_1_M .* (x1(:, p,:)-x_mrkwide(:,p,:)));
    expanded_sum_product_1_M = repmat(sum_product_1_M, J, 1, 1);

    sum_product_0_M = sum(share_0_M .* (x0(:, p,:)-x_mrkwide(:,p,:)));
    expanded_sum_product_0_M = repmat(sum_product_0_M, J, 1, 1);

    GradDelta_Y_M(:,j,:,k,p) = share_1_M.*(x1(:,p,:)-x_mrkwide(:,p,:)-expanded_sum_product_1_M)-share_0_M.*(x0(:,p,:)-x_mrkwide(:,p,:)-expanded_sum_product_0_M);
end
if model == "RCL"
for i = 2:3
    dtheta_1 = f_calc_gradient_mktwide( ...
        xnonlin_mrkwide(:,i,:), ...
        dfull_sim_mrkwide,   ...
        share_main_1_ijt_M, ...
        dDelta_RCL(:,i,:) ...
    );

    dtheta_0 = f_calc_gradient_mktwide( ...
        xnonlin_mrkwide(:,i,:), ...
        dfull_sim_mrkwide,   ...
        share_main_0_ijt_M, ...
        dDelta_RCL(:,i,:) ...
    );

    GradDelta_Y_M(:,j,:,k, length(theta_M)+i-1) = dtheta_1 - dtheta_0;
end
end

end
end

% agg to brand level
GradDelta_Y_M_agg   = zeros(J,J,nmarket,K,size(param_M,1));
GradDelta_Y_M_agg(:,:,:,:,1) = GradDelta_Y_M(:,:,:,:,1);
for i = 2:14
    GradDelta_Y_M_agg(:,:,:,:,i) = sum(GradDelta_Y_M(:,:,:,:,3*(i-2)+2 : 3*(i-2)+4),5);
end
if model == "L_brandfe_montfe" || model == "RCL"
    GradDelta_Y_M_agg(:,:,:,:,15:end) = GradDelta_Y_M(:,:,:,:,41:end);
end

F = 1/(nmarket*K)*sum(sum(sum(sum((Delta_Y-Delta_Y_M).^2))));
dF = zeros(length(param_M),1);
for p = 1:size(param_M,1)
dF(p) = -1/(nmarket*K)*sum(sum(sum(sum(2.*(Delta_Y-Delta_Y_M).*GradDelta_Y_M_agg(:,:,:,:,p)))));
end

end